%{
AUTHOR: Felipe Arteaga
-------------------------------------------------------------------------
PROJECT: Warnings
-------------------------------------------------------------------------
DESCRIPTION:
=========================================================================
%}


clearvars -except projectDir projectDirData fromMainWarningsPaper
clc;close all;fclose('all');feature('DefaultCharacterSet','UTF-8');

if(not(exist('projectDir','var')==1&&exist('projectDirData','var')==1&&exist('fromMainWarningsPaper','var')==1&&fromMainWarningsPaper))
    pcName=char(java.lang.System.getProperty('user.name'));
    if(strcmp(pcName,'felipe'))
        % PC Felipe
        myDir='/Users/felipe/Dropbox/';
        projectDir=[myDir,'git/warnings/'];
        projectDirData=[myDir,'projects/warnings/'];
        addpath(genpath([myDir,'/myMatlabFunctions/']));
    end
end
compileLatexTable=false;

%%

dirPlots=[projectDir,'/paper/figuresCL/RDs/'];

% warning('Alternative dir plots')
% dirPlots=[projectDir,'/paper/figuresCL/RDs/Aux/'];

dirPlotsR=[projectDir,'/paper/figuresCL/RDs_multBandwidth/'];
dirTable=[projectDir,'/paper/tablesCL/'];

dirData=[projectDirData,'/data/chile/'];


closePlotsWhilePlotting=true;

plotRDs=false;
savePlotRDs=true;
plotWithZoom=true;
bandwidthTable='local'; % 'full','local','localAuto';


% For making subfigures with plots:
maxHorzPlots=3;
maxVertPlots=4;

tableBandwidthComparison=false; % Make a table comparing estimates with different bandwidths.
posBandwidth=1; % This is for the bandwidht robust excercise. Indicates the subpopulation that we are using the bandwidth table
plotAndSaveRobust=true; % Plots of robust bandwidht RD

posDetails=1;% This is for the details RD table. Indicates the subpopulation that we are using for the details table

switch bandwidthTable
    case 'full'
        bwnote='Full bandwidth was used (-0.28,+0.68) with 2nd order (local) polynomial fit. ';
    case 'local'
        bwnote='Bandwidth was set to +-0.1 with default local polynomial fit (1st order). ';
    case 'localAuto'
        bwnote='Bandwidth was set to automatic (bwselect with default options) with default polynomial fit (1st order). ';
end

%% Load data:
anho=0; %0: pooled
if(anho==0)
    anhoStr='';
else
    anhoStr=sprintf('%i',anho);
end


notIn2020={'enrolledInAssigned','enrolledInAddedIniEnd','enrolledInAddedNoWaitlistIniEnd','enrolledInAddedNoRiskIniEnd'}; % {'assignedToPref','acceptOffer','declineOffer','preferenciaAsign_1','placedInAnyPreferenceIniEnd','placedInAddedIniEnd','placedInAddedNoWaitlistIniEnd','placedInOldIniEnd','enrolledInAssigned','enrolledInAddedIniEnd','enrolledInAddedNoWaitlistIniEnd'};
if(anho==0)

    d1=load([dirData,'2018/inputRD'],'dataRD');
    d2=load([dirData,'2019/inputRD'],'dataRD');
    d3=load([dirData,'2020/inputRD'],'dataRD');
    % Variables that are not available for 2020 yet:
    fillWithNan=notIn2020;
    %fillWithNan=
    for v=1:length(fillWithNan)
        d3.dataRD.(fillWithNan{v})=nan(height(d3.dataRD),1);
    end


    % Universe that is comparable between 2018-2020
    entryGrades=[-1 0 1 7 9];
    sirven1=ismember(d1.dataRD.grade,entryGrades)&not(d1.dataRD.cod_reg==13);
    sirven2=ismember(d2.dataRD.grade,entryGrades)&not(d2.dataRD.cod_reg==13);
    sirven3=ismember(d3.dataRD.grade,entryGrades)&not(d3.dataRD.cod_reg==13);

    d1.dataRD.comparable=sirven1;
    d2.dataRD.comparable=sirven2;
    d3.dataRD.comparable=sirven3;

    d1.dataRD.anho=2018*ones(height(d1.dataRD),1);
    d2.dataRD.anho=2019*ones(height(d2.dataRD),1);
    d3.dataRD.anho=2020*ones(height(d3.dataRD),1);

    % Keep vars in common

    incommon=intersect(intersect(d1.dataRD.Properties.VariableNames,d2.dataRD.Properties.VariableNames),d3.dataRD.Properties.VariableNames);
    incommon=incommon(not(strcmp(incommon,'mrun')));
    dataRD=[d1.dataRD(:,incommon);d2.dataRD(:,incommon);d3.dataRD(:,incommon)];



else
    load([dirData,anhoStr,'/inputRD'],'dataRD')
end

dataRD.rural=dataRD.newMarket<0;

% As asignment at .3 is control, just move (if any) in .3 to epsilon to the
% left:
dataRD.riskPopup(dataRD.riskPopup==.3)=.29999999;

% Avoid mass points at extremes:
restrAll=dataRD.pobPopup==1&dataRD.riskPopup>.01&dataRD.riskPopup<.99;

% No need to be more specific, and put "Predicted placement risk of 1st
% attempt"
dataRD.Properties.VariableDescriptions{'riskPopup'}='Predicted placement risk';
dataRD.Properties.VariableDescriptions{'treatedPopup'}='Treated pop-up in first attempt';


%%
qsCell=cell(3,3);

qsCell{1,1}=discretize(dataRD.totalEEsIn2kmR_geo,[0,quantile(dataRD.totalEEsIn2kmR_geo,.1:.1:.9),105],'categorical');
qsCell{1,2}='Schools in 2km radio from geocoding';
qsCell{1,3}='totalSchoolsIn2km';

qsCell{2,1}=discretize(dataRD.totalEEsIn3kmR_geo,[0,quantile(dataRD.totalEEsIn3kmR_geo,.1:.1:.9),105],'categorical');
qsCell{2,2}='Schools in 3km radio from geocoding';
qsCell{2,3}='totalSchoolsIn3km';

dataRD.cantEECat=renamecats(dataRD.cantEECat,{'Santiago (>700)'},{'Santiago ($>$700)'});

qsCell{3,1}=dataRD.cantEECat;
qsCell{3,2}='Market size';
qsCell{3,3}='marketSize';


for k=1:size(qsCell,1)
    for o=1:3
        switch o
            case 1

                outcome='addAnyIniEnd';
                labelOutcome='Add any school';
            case 2
                outcome='DRiskExpostIniEnd';
                labelOutcome='$\Delta$ Risk';
            case 3
                outcome='schoolsAddedIniEnd';
                labelOutcome='Schools Added';
        end
        qs=qsCell{k,1};

        uniqueQs=unique(qs);
        uniqueQs=uniqueQs(not(ismissing(uniqueQs)));
        cantQs=length(uniqueQs);


        matrix_b=nan(cantQs,1);
        matrix_sd=nan(cantQs,1);
        for q=1:cantQs

            pop=restrAll&qs==uniqueQs(q);
            res=rdrobust(outcome,'riskPopup',dataRD(pop,:),'subpop',true(sum(pop),1),'cutoff',.3,'rdrobustOpts','h(.1 .1) p(1)','withPlot',false);
            %title(sprintf('nm: %i , year: %i',nm,a))
            % Beta
            matrix_b(q)=res.tau_cl;
            % Sd Beta
            matrix_sd(q)=res.se_tau_cl;
        end
        figure
        bar(uniqueQs,matrix_b)
        hold on
        errorbar(uniqueQs,matrix_b,1.96*matrix_sd,'lineStyle','none')
        xlabel(qsCell{k,2});
        ylabel(sprintf('RD estimate - %s',labelOutcome))
        easyExport([dirPlots,'rdEstimates_',outcome,'_',qsCell{k,3}])

    end
end
